Instability of Boost-invariant hydrodynamics with a QCD inspired bulk viscosity 
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We solve the relativist ic Navier-Stokes equations with homogeneous boost- invariant boundary 
conditions, and perform a stability analysis of the solution. We show that, if the bulk viscosity has 
a peak around Tc as inferred from QCD-based arguments, the background solution "freezes" at Tc 
to a nearly constant temperature state. This state is however highly unstable with respect to certain 
inhomogeneous modes. Calculations show that these modes have enough time to blow up and tear 
the system into droplets. We conjecture that this is how freeze-out occurs in the QGP created in 
heavy ion collisions, and perhaps similar transitions in the early universe. 
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PACS numbers: 25.75.-q,25.75.Dw,25.75.Nq 

The system produced at RHIC []?] is believed to be a 
good liquid in the early stages of its evolution. At a later 
stage this liquid transforms into a gas of particles which 
interact more weakly and eventually decouple. 

The transition from a strongly interacting liquid into 
particles is however still not understood, both on the fun- 
damental and phenomenological level. On a conceptual 
level, no adequate description exists of how the mean- free 
path goes from zero ( "nearly ideal liquid" ) to a distance 
comparable to system size ("transport" of particles). 

On a phenomenological level, hydrodynamics fails to 
reproduce particle interferometry data [2] unless the sys- 
tem decouples at an unrealistically high temperature ^] 
Tc, the critical temperature for the QCD phase tran- 
sition) . Moreover, attempts to use a hadronic transport 
model as an afterburner to hydrodynamics (conceptu- 
ally considered to be the "next best" approximation) has 
failed to improve the model-data agreement [H, 

Recent attempts to reconcile hydrodynamics with in- 
terferometric data have focused on a supposed sharp 
increase, and perhaps divergence, of bulk viscosity near 
Tc. This behavior of bulk viscosity has now been inferred 
from a variety of arguments [lOl, lUj. In this work, we 
combine the recently acquired understanding of viscosity 
with an earlier study [12] to get a simple picture of how 
the peak in viscosity triggers freeze-out process. 

The Navier-Stokes equations with Boost-invariant 
symmetry [l3|, [l5| can be rewritten [T2[ in terms of the 
Reynolds number i?, the entropy s the co-moving time 
r, the total number [N) of dimensions, and the dimen- 
sionality of the homogeneous expansion (M) 
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For example, M = 1 = 3 corresponds to the case 
studied in [13], M = AT = 3 to the "Krakow model" [3 
or to a Friedman-like solution in flat space. 

The Reynolds number is a function of temperature 
T,bulk and shear viscosity and r] and entropy s 
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2(1 - M/N)r] + MC 

TST 
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Given expressions for 5, 77, in terms of T, this set of 
equations becomes closed and solvable. 

For the equation of state, we use the parameterization 
of the speed of sound given in [16] (We have checked 
that our results do not vary qualitatively if the ideal EoS 
is used). We use the "sCFT limit" [I^ for r] {r] = s/Aii) 
and a parametrization of [TO] for 



C = -5 ZpQCD 
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where t = T-T^ and a = O.OlTc and z^qcd ^ 10"^ [ij. 
At T > Tc (through not at the hadronic regime, T < Tc)^ 
this ansatz provides a reasonable fit to [9]. 

We follow the stability analysis performed in [l2[ • The 
amplitude of a generic perturbation to the ID Boost- 
invariant system is a vector x in the two dimensional 
space of entropy perturbations and flow (rapidity y) per- 
turbations 
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and its dependence on rapidity can be decomposed into 
Fourier components in rapidity space of wavenumber k 



ky 



(5) 



The equation of motion for x will then be given by 
d f xi \ f All A12 \ f xi 



dr V 
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where Aij is a real matrix function of and k. We 
refer the reader to Eq. 4.23-4.26 of ^ for the full form 

of Aij . The evolution of perturbations can be understood 
through the behavior of the modulus of X = x^x 



T—X x^ MijX 
dr 



(7) 
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Since Mij is real and symmetric, it will always have two 
real Eigenvalues, Xmax and Xmin (corresponding Eigen- 
vectors Xmax,min)i as wcll as orthogoual matrices Bij di- 
agonalizing it. Defining yi = B^j Xjwe see that 

2 2 /r>\ 

Xminymin ^ ^ XmaxUmax \^) 

Thus, if Xmin > 0, the system is inherently unstable, 
since perturbation in any direction will produce a positive 
growth rate. If Xmax < 0, on the other hand, the system 
will be stable against all perturbations. 

If Xmin < and Xmax > somc modes will be stable 
and some will be unstable. In the latter case, the time de- 
pendence of Aij will in general continuously rotate Xmin 
and Xmax in time: While an initial perturbation in the 
Xmax direction grows as a power-law with r/rini^ where 
Tini is the starting time of the perturbation 

X{r) ^ ( — j Xmax{rini) (9) 

\ '^ini / 

its growth might well be stopped by a change of direction 
of the eigenvalues Aij{r): Under a general evolution of 
Aij J Xmaxi^ini) might be in the direction of Xmin 

(r > 

Tini) a short time later. Solving Eq. [6] will take this 
effect into account. 

In what follows, we use zq = 0,0.1,1 Tc and evolve 
the system from an initial temperature T = 0.3 GeV 
and comoving time 0.6 fm. Note that the qualitative 
features of our study are independent of the details of 
the evolution before Tc, such as the initial temperature 
and timescale. 

Fig. [1] shows the temperature, entropy density and 
total entropy in the central rapidity unit as a function 
of time. As can be seen, as soon as zq becomes non- 
negligible (i.e. viscous forces dominate around Tc), the 
kinematic evolution of the system "freezes" . The system 
then stays at nearly constant temperature, through it 
keeps producing more and more entropy at the expense 
of advective energy. 

At first sight, large values of zq are excluded by HBT 
data and multiplicity measurements However, we 

will show that this long phase is unstable against small 
perturbations. Thus, its further evolution will not be 
given by the background solution, but by a rapid forma- 
tion of local inhomogeneities. 

Fig. [2] shows the Xmin and Xmax Eigenvalues corre- 
sponding to representative /c = 2, 8 (other values of k 
were checked not to vary significantly wrt those presented 
here). As can be seen, the peak in bulk viscosity forces 
the growing/damping rates to increase rapidly. Thus, 
any initial perturbation in the unstable direction will 
rapidly grow to a value comparable with the background, 
unless the system's evolution will stop the growth by ro- 
tating the direction of the unstable modes. 



Fig. [3] examines weather this occurs for larger values 
of Zq. If the peak of viscosity is negligible, the unstable 
eigenvector keeps rotating throughout the evolution of 
the system. Thus, even an unstable mode's growth will 
very quickly stop growing since the dynamics will turn 
it into a damped mode. When zq dominates, however, 
something every interesting happens: At the time when 
T approaches Tc, the direction of the unstable modes 
experiences an abrupt rotation. Then it stays constant 
throughout the time the system travels through the vis- 
cosity peak (this time increases strongly as zq increases), 
and gets rotated again as the peak is passed. The rea- 
son for this behavior is clear: Fig [2] shows that the peaks 
occur at T Tc, where the background of the system is 
dynamically "frozen" (all advective energy is turned into 
entropy by viscous processes to keep temperature nearly 
constant). Thus, it is not surprising that the direction 
of the unstable Eigenvector also remains approximately 
constant. 

Thus, at large zq unstable modes have all the required 
time to grow, especially considering the growth rate, ~ 
Xmax in Fig. [21 becomes overwhelming compared to the 
background expansion rate. 

Finally, Fig. [H shows the explicit solution of Eq. [6l 
At each time-step, a perturbation is born in the unstable 
Eigenvector mode, and then evolved until the end of the 
evolution of the system. The plot shows X (r) / X (rini) ^ 
the ratio of strength of the perturbation to the initial 
strength as a function of time (Note that the starting 
value is always unity). A large X{r)/X{rini) does not 
mean the evolution equation is invalid: For any point in 
the graph, there will be a small enough perturbation am- 
plitude that survives as a perturbation in the subsequent 
evolution. The probability of a larger perturbation form- 
ing and significantly modifying the background, however, 
should grow strongly with zq. 

Fig. [4] makes it clear that any microscopic mechanism 
seeding instabilities at the scale X ~ 10~^,uniformly dis- 
tributed in X and at a rate of ~ fm~^ is likely to generate 
power-law growing instabilities a few fm after T ^ Tc. 
These instabilities should reach X > 1, and hence play an 
important role in the subsequent evolution of the system, 
a few fm after that. 

Our model does not have the required scope to tell 
us what this role is likely to be. Once grown, perturba- 
tions will break local homogeneity and Boost-invariance, 
requiring full 3D hydro to be modeled properly. 

Transverse flow also has the potential to modify this 
scenario. We note, however, that substituting ID for ho- 
mogeneous 3D Boost-invariant hydrodynamics [12] does 
not quantitatively change our results. As this "3D 
Bjorken equation" (corresponding to the hydro-inspired 
model examined in [ij]) is locally similar (up to rescaling) 
to the asymptotic behavior of realistic solutions describ- 
ing transverse expansion [205, trust that the effect 
pointed out in this work is relevant for systems without 
transverse homogeneity. 

While we leave the stability analysis in second order 
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FIG. 1: (Color online) Thermal parameters as a function of time for zo — 0,0.1, ITc Solid line shows T /Tc. Dot-dashed shows 
entropy density dS/dy normalized to initial value, and Dashed the integrated entropy S from an initial unit of rapidity, again 
normalized to the initial value. 
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FIG. 2: (Color online) Xmin and Xmax as a function of T /Tc^ for zq = 0,0.1, ITc. Only mixed stability data-points are shown 
(at later times in the = case, evolution becomes stable). Solid black corresponds to = 2,dashed blue to A; = 8 



hydrodynamics (Israel- Stewart [2ll , [22I , [23| or its modern 
variations [H, [H, [13, S (3) to a future project, 
we do not expect it to play a big role in starting the 
instabilities: As argued in [8], as long as the system's 
expansion rate as the system approaches Tc is smaller 
than the relaxation time rn 

m'-f<^ (10) 

second order terms will not prevent viscous corrections to 
the pressure of the order of C/t, but merely localize their 
propagation. The (admittedly unreliable) estimates from 
strongly coupled CFT [30] suggest that this criterion is 
amply satisfied, especially considering that the increase 
in viscosity causes the background solution to slow down 
over a timescale much bigger than rn- 2nd order hydro- 
dynamics might however play a dominant role in stopping 
the instabilities once they grow to a value comparable to 
the background. 

Viscosity and causality should prevent any subsequent 
reinteraction of the created inhomogeneities. It is highly 
likely, therefore, that instabilities evolve isolated frag- 
ments, moving away from each other with pre-existing 
longitudinal and transverse flow. This will have the ef- 
fect of quickly stopping the unrealistically high entropy 
production seen in Fig. [H 



The small dependence on k of the instabilities growth 
rate (Fig. [4]) might cause the formed inhomogeneities to 
continue seeding smaller and smaller instabilities, until 
the size of the smallest inhomogeneity becomes ~ ^qcd^ 
the scale at which local expansion will be slowed down 
by the viscous forces. The temperature of the center of 
such a "cluster" should be ~ Tc. It is therefore tempting 
to regard these clusters, once formed and decoupled, as 
the "fireballs" within the Hagedorn picture [31]. 

This scenario shares some phenomenological similar- 
ities with [3I, m, [35I while differing in the funda- 
mental description: In [32|, [sl, [35[,The role of the 
instability generator is the first order coexistence phase, 
and the instabilities are defined by being in a phase differ- 
ent from the "background" . In the current scenario, the 
instabilities are generated purely through hydrodynamic 
evolution, no discontinuities in the Equation of state are 
needed, and the perturbations are purely hydrodynamic 
in origin. 

While this work focuses on heavy ion collisions, the dy- 
namics examined here could play a role in cosmology, pro- 
vided that quantum fields involved in cosmological phase 
transitions are approximately thermalized, and their vis- 
cosity depends on T in a similar way to the one assumed 
here. The general relativistic solution could than be lo- 
cally similar enough to the (3D) Boost-invariant expan- 




sion (as in [36]) to allow for formation of local instabili- 
ties. These could contribute to structure formation, pro- 
vide the loci where baryogenesis would occur |37| , or seed 
the formation of microscopic black holes. The non-linear 
attractive behavior of gravity could enhance instability 
growth above hydrodynamic expectation. 

In conclusion, we have shown that, when the conjec- 
tured rise of bulk viscosity at Tc is inserted into the 
Navier- Stokes equations with Boost-invariant symmetry, 
the solution "freezes": While entropy continues to in- 
crease, temperature remains constant. The background, 
however, becomes unstable against small perturbations. 



which then grow until boost-invariant hydrodynamics 
breaks down . We hope that further work will clarify 
the subsequent evolution of these instabilities, and pro- 
vide a quantitative link between this scenario and data. 
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